Ordering and Excitation in Orbital Compass Model on a Checkerboard Lattice 
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We study an orbital compass model on a checkerboard lattice where orbital degree of freedom is represented 
by the pseudo-spin operator. Competition arises from an Ising interaction for the z component of pseudo- spins 
along the vertical/horizontal bonds and an Ising interaction for the x component along diagonal bonds. Classical 
and quantum compass models are analyzed by utilizing several analytical methods and numerical simulations. 
At a fully frustrated point where the two Ising interactions compete with each other, a macroscopic number 
of orbital configurations are degenerate in a classical ground state. This degeneracy is lifted by thermal and 
quantum fluctuations, and a staggered long-range order of the z component of the pseudo-spin is realized. A 
tricritical point for this order appears due to competition between the bond dependent Ising interactions. Roles 
of geometrical frustration on excitation dynamics are also examined. 

PACS numbers: 75.25.Dk, 75.30.Et,75.47.Lx 
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I. INTRODUCTION 

Long-range order and excitation dynamics in the orbital de- 
generate correlated electron systems are one of the recent at- 
tractive themes in condensed matter physics.^ Orbital degree 
of freedom represents a spatially anisotropy of the electronic 
wave function. In molecules, orbital degeneracy is usually 
lifted by coupling with lattice, i.e. the Jahn-Teller effect. In 
contrast, in crystal lattices, there are some equivalent bonds 
around a transition-metal ion. When an orbital is directed 
along one of the equivalent bonds, anisotropy in the bond 
energies comes out. In this sense, all bond energies on the 
equivalent bonds are not minimized simultaneously. This is 
regarded as a kind of frustration effect termed "orbital frustra- 
tion"i2i^ This characteristic in the orbital degenerate systems 
provides a wide variety of exotic phenomena such as order by 
disorder phenomena^^ orbital liquid state^'^ and so on. 

One of the well- studied orbital models is the Kugel- 
Khomskii model. ^ This is applied to orbitally degenerate Mott 
insulators, where the orbital degree of freedom in a transition- 
metal ion is described by the pseudo-spin (PS) operator. The 
intersite interactions between the nearest neighbor ions are 
represented by products of the Heisenberg-type interaction be- 
tween spins and the orbital interaction. In the orbital part, 
the interaction between the pseudo-spins explicitly depends 
on the bond direction. 

Another well studied orbital model is the orbital compass 
model where the orbital degree of freedom is only taken into 
account. A general expression of the orbital compass model 
is given by 



H 






Ti){ni-Tj), 



(1) 



where T^ is the PS operator for the doubly-degenerate orbital 
degree of freedom with an amplitude of 1/2, fii is a unit vec- 
tor along the bond direction /, and {ij)i indicates the nearest 
neighboring (NN) i and j sites along /. This model has some 
analogy to the dipole-dipole interaction and shows the same 



characteristics with the orbital part in the Kugel-Kohmskii 
model; the interactions explicitly depend on a bond direc- 
tion. This model has been studied from broad view points; 
quantum phase transitions^ topological quantum order^\ hid- 
den dimer order,i^ and protected qubit^^'^^ are examined on a 
square-lattice compass model, and a kind of compass model is 
proposed as an appropriate model for transition-metal oxides 
with a strong spin-orbit coupling. ^^ 

Transition-metal ions with orbital degree of freedom some- 
time consist geometrically frustrated lattices, such as trian- 
gle and spinel crystals. Interplay of geometrical frustration 
and spin-orbital entanglement often give rise to novel states of 
matter, such as spin-orbital molecules in AlV204r^ a cooper- 
ative release of frustration proposed in ZnV204r^ and a res- 
onating valence bond state predicted in LiNi02 J^ Even with- 
out spin-orbital entanglement, ground state and excitation dy- 
namics in orbital degenerate system with geometrically frus- 
trated lattice are non-trivial because of the orbital anisotropy 
and frustration characteristics. The recent neutron scattering 
experiments suggest an excitation from spin-orbital molecules 
in GeCo204,^^ where effective total angular moments might 
be described by the orbital compass model. 

In this paper, the orbital compass model on one of the ge- 
ometrically frustrated lattices, i.e. a checkerboard lattice, is 
studied. Competitions arise from an Ising interaction for the 
z component of PS, T^, along the horizontal and vertical di- 
rections on a lattice and an Ising interaction for T^ along the 
diagonal directions. Phase diagrams in classical and quantum 
models, where PS operators are regarded as classical vectors 
and quantum operators, respectively, are obtained by several 
analytical and numerical methods. It is shown that, at a fully 
frustrated point, where a number of classical PS configura- 
tions are degenerate, two-dimensional staggered T^ ordered 
state is stabilized by thermal and quantum fluctuations. Be- 
cause of the bond depend Ising interactions, a tricritical point 
for the two-dimensional staggered T^ order appears. A one- 
dimensional characteristic excitation in this ordered state is 
remarkable near the phase boundary. Present results are com- 
pared with the results in the compass model on a square lat- 



tice. 

In Sect. II, an orbital compass model on a checkerboard lat- 
tice is introduced. In Sect. Ill, a classical model is analyzed by 
the mean-field (MF) approximation and the classical Monte- 
Carlo (MC) simulation. In Sect. IV, a quantum model is an- 
alyzed by the spin-wave approximation, a combined method 
of the MF approximation and the Jordan- Wigner transforma- 
tion, the exact-diagonalization method and the quantum MC 
simulation. In Sect. V, results for the excitation dynamics are 
presented. Section VI is devoted to summary and discussion. 



II. MODEL 

We set up an orbital-compass model on a two-dimensional 
checkerboard lattice. Doubly-degenerate orbitals, dzx{= o) 
and dyz{= b), are introduced in each site. We focus on the 
orbital degree of freedom of electrons and neglect the spin 
degree of freedom. We start from the spinless Hubbard model 
with the doubly degenerate orbitals defined by 



<ij>i 



<ij>'rn 



^E 



'^ia'^ihi 



(2) 



where c^cr is the annihilation operator for a spinless fermion 
with orbital 7(= a, 6) at site z, and symbols < ij >i and 
< ij >^ represent the NN- and the next NN (NNN) ij pairs, 

respectively. The transfer integral t^^^.^^, ('^n'nn^jY^ is de- 
fined on a NN (NNN) bond along the direction / (m). Ma- 
trix elements of the transfer integrals are determined by the 
Slater- Koster parameters. Since an electron in the dzx (dyz) 
orbital hops along the x {y) direction, the matrix elements for 
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with a positive constant ti, where a are the Pauli matrices. 



Matrix elements for t 
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linear combinations of the dz 



dyz)/V2, 



/ are obtained by introducing the 
and dyz orbitals, i.e. {dzx i 



as 
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with a positive constant ^2 . From this Hubbard-type Hamil- 
tonian, an effective Hamiltonian in the case ofU ^ ^1,^2 is 
derived by the second-order perturbational procedure as 



H = J, 
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(7) 



where we define the exchange constants Jz = 2t1/U and 
Jx = 211/11 . We introduce the PS operator, T^, with a magni- 
tude of 1/2, where the dzx and dyz orbitals are taken to be the 
eigen states of Tf. This Hamiltonian is a kind of the orbital 
compass model defined on a checkerboard lattice in a sense 
that the Ising-type interactions depend on bond directions. 

Next we discuss a symmetry of the Hamiltonian in Eq. d?]). 
Let us focus on a NNN bond network on a checkerboard lat- 
tice (see Fig. [T]). One dimensional chains along (11) and 
(11) directions are independent with each other. In one of 
the chains, termed /, we introduce the operator defined byi^ 



iei 



(8) 



where i runs along this chain. It is shown that this operator 
commutes with the Hamiltonian by using the commutation 
relation [a^a^^afaj] = with i ^ j. Therefore, the en- 
ergy eigenstates are labeled by the eigenvalues of Pi, i.e. ±1. 
There are L labels on a L x L-site lattice. This character- 
istic is available in numerical exact-diagonalization calcula- 
tions for large cluster size. Because of these local symme- 
tries, the generalized Elitzur's theorem is applicable to this 
model*^ It is rigorously shown that a long range order of 
Tq = N~'^ ^^ j^x^iq-vi ^Qj. ^^y momenta of q, which does 
not commute with Pi, is not realized at finite temperature. 



III. CLASSICAL ORBITAL STATE 



FIG. 1: A schematic picture for the orbital compass model on a 
checkerboard lattice. 



In this section, we treat the orbital PS as a classical vector 
defined in a two-dimensional T^-T^ plane with an amplitude 
of 1/2. 
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FIG. 2: Phase diagram obtained by the MF approximation. The 
Unes in J < 1 and J > 1 are plotted in the different scales at left 
and right figure, respectively. Stable orders are (T^^^^^) for J < 1, 
and (T(^ Q)) or (T^q^^) for J > 1. Transition temperatures do not 
depend on J and are Tc/{2Jz) = 0.5 in J < 1 and Tc/Jx = 0.5 
in J > 1. Insets show schematic PS configurations in the T^ — T^ 
plane. 
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FIG. 3: Phase diagram obtained by the classical MC method. The 
data in J < 1 and J > 1 are plotted in the different scales at left and 
right figure, respectively. Bold and double lines represent the second- 
and first-order phase transitions, respectively. Dotted line represents 
the crossover below which an one-dimensional T^ correlation devel- 
ops. The filled circle at J = oo is determined by applying the MC 
simulation to the one-dimensional Ising model on L = 200 chain. 
The first order occurs at J = 1 and zero temperature. 



A. Mean-Field Analysis 

First, we show the orbital state obtained by the MF approx- 
imation. We take the MFs for the orbital order as (T^) = 
AT-i Ei(^Oe'^''^ for (/ = x,z). In Fig. [2 the phase dia- 
gram in the plane of J = Jx/{2Jz) and temperature, T, is 
presented. Stable orbital orders are (Tf ^x) for J < 1, and 
(Tf Q^) or {T?Q ^^) for J > 1, i.e. the staggered T^ order 
along (11) and (11) directions. The transition temperatures 
do not depend on a magnitude of J. Beyond the analyses 
for the MF order parameters with single momentum, there 
are a number of degenerate MF solutions for J > 1 ; in the 
(T.Q ^ ) ordered state, we consider the transformation of PS 
that T^ -^ —T^ for all sites in a certain chain along (11) 
and (11) directions. The MF energy is not changed under this 
transformation, since T^ operator is only concerned in the in- 
teraction along the (10) and (01) directions. There are 2^^ 
degenerate MF solutions on a L x L-site lattice at T = 0. 

At a point of J = 1 and T = 0, there is an additional 
degeneracy. Any linear combinations of (Ti^ ^^) and (Tf ^^), 
i.e. {T{0)) = cosO{T^l^^^) + sin^(T(-o^^)') where is the 
rotation angle in the T^ — T^ plane, have the same energy. 
This degeneracy is not expected from the Hamiltonian which 
does not show any continuous symmetry. 



B. Monte Carlo Simulation 

In this subsection, we introduce the numerical results ob- 
tained by the classical MC simulations. Two-dimensional 
20^-, 30^- and 40^-site clusters with a periodic boundary con- 
dition are used. We adopt the Wang-Landau algorithm^ 



where 5 x lO^MC steps are used for both making histograms 
and measurements. 

The phase diagram is presented in Fig. [3] With increasing 
J from J = 0, where the model is reduced to the T^-Ising 
model on a square-lattice, the transition temperature for the 
(Tf ^ ) order gradually decreases because of the competition 
between J^ and Jx. At another limit, J = cx), the model is 
reduced to the independent one-dimensional T^ -Ising model 
which does not show a long-range order. However, there is a 
crossover temperature around T/Jx = 0.1 where the specific 
heat C shows a broad peak below which an one-dimensional 
T^ correlation develops. By introducing the NN interaction, 
J^, the broad peaks remain in the temperature dependences 
of the specific heat. The peak positions are plotted by dot- 
ted lines in Fig. [3l The crossover temperature gradually de- 
creases, when the system approaches to the J = 1 point. Let 
us focus on the point of J = 1. At T = 0, continuous de- 
generacy exists as explained above. With increasing T, the 
(Tf ^^ ) order is stabilized among them due to the thermal ef- 
fect. This is a kind of order by fluctuation phenomena. When 
temperature is increased furthermore, a conventional thermal 
effect makes a system to be an orbital disordered state. As 
a result, reentrant feature is observed in the phase boundary 
around J = 1. 

Next, we present detailed MC results around J = 1. 
In Fig. m the temperature dependences of the specific heat 
and the squre root of the staggered PS correlation function 
S'^^(7r, tt) are presented on several cluster sizes at J = 1. We 
define S''{q) = N'^ ^.. i^^T/e^^^^^-^^). A sharp peak is 
observed in C at T/(2Ja;) ~ 0.12 which is termed Tc from 
now on. A value of C = 0.5 in the limit of T = implies 
an existence of one degree of freedom per site, i.e. the po- 
lar angle of PS in the T^ — T^ plane. The correlation func- 




FIG. 4: Temperature dependences of (a) specific heat and (b) square 
root of the staggered-type correlation function 5'^^(7r, tt) for several 
cluster sizes at J = 1. 



tion starts to increase around Tc and approaches to the upper 
limit of 0.5 at low temperatures. With increasing the system 
size, Tc slightly decreases, a peak in C becomes sharp, and 
an increase in 6'^^(7r, tt) at Tc becomes sharp. The results im- 
ply that these anomalies at Tc correspond to the second-order 
phase transition in the thermodynamic limit. 

In Fig. \5\ the numerical results at J = 1.045 are pre- 
sented. We show the temperature dependences of energy E, 
C, [5'^^(7r,7r)]i/2 ^^^ {T^f^'^ (I = x,z) for several size 
clusters. Two anomalies are observed diiT/(2Jz) ~ 0.067 
and 0.09 which are termed Tl and Th, respectively. At Th, 
a peak in C becomes sharp, and an increase in 5'^^(7r, tt) be- 
comes remarkable with increasing the system size. These re- 
sults are similar to the results at Tc in J = 1 (see Fig. [4]). 
On the other hand, at Tl, S^^{7r, tt) and (T^^) decrease and 
(T^^) increases with decreasing T. Large system- size depen- 
dences are observed in E and C at T^. These results imply 
that, in the region of Tl <T < Th, the {T^^^^^) order is real- 
ized in the thermodynamic limit. Below Tl, PS's are directed 
along the T^ axis. This is expected from the one-dimensional 
T^ Ising interaction. 

To examine the orders of the phase transitions at Tl, we 
calculate the free energy defined by 



F = -T\n"^D{E)e 
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FIG. 5: Temperature dependences of (a) energy, (b) specific 
heat, and (c) square root of the staggered-type correlation function 
S^^{q = tt), root mean squares ofthe PS moment (7^^^), and (IT^) 
for several cluster sizes at J = 1.045. 



where D{E) is the density of states. We calculate y{E) = 
\nD{E) — (3 EN -\- const, as the energy histogram in the 
Wang-Landau scheme in the MC simulation.^^ The results for 
three temperatures around Tl are shown in Fig. [6] as func- 
tions of E. Double peak structures are commonly observed 
in y{E). The two peak heights are reversed by changing tem- 
perature; y{E) at the lower-energy peak increases with de- 
creasing T. Since y{E) is proportional to the system size 
as expected from the definition, the two minima in —y{E) 
are separated by energy of the order of TV. It is expected in 
the thermodynamic limit that the energy in the stable state is 
changed discontinuously by changing T, and the anomaly at 
Tl corresponds to the first-order phase transition. 



IV. QUANTUM ORBITAL STATE 

In this section, the orbital PS's are treated as quantum spin 
operators with an magnitude of 1/2. 
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FIG. 6: Energy dependence of ^(£;) = \\iD{E) - /3EN + const, 
(see text) for several temperatures at J = 1.045. The cluster size is 
chosen to be A^ = 40^. 



A. Spin- Wave Approximation 

The ground state at T = is analyzed by using the 
spin wave approximation. The long-range ordered states of 
{T? ^.) and (T.^^^) are adopted as the ground states for 
J < 1 and J > 1, respectively, although a number of de- 
generate states exist for J > 1. At J = 1, continuous PS 
configurations connecting (T^^^^)) to (T^^^,,)), i.e. {T{lp)) = 
cos ^{T? ^) + sin (p{T?Q ,) are assumed (see Fig. [8]). There 
are four sublattices in the ordered states. 

By introducing the four kinds of the Holstein-Primakoff 
bosons, (ttfc, bk^Ck, dk), the Hamiltonian up to the second or- 
der of the boson operators is given as 

Hsw/{2Jz) = -NS^icos^ LP + Jsin^ Lp) 

N/A 

+ 5'^ 2 (cos^ LP + Jsin^ Lp) hi 



k >- 
-hlsin^ip^-hlcos^if 



(10) 



where S = 1/2, and 

^k = (4^fc + ^k^k + clck + dldk), 

hi = cos kx {aibk + cldk + ^l^Lfc 4 

- cos ky {a'lck + b]^dk + ct'lcij^ -] 



(11) 



■cidik- 

-bidlk 



H.c.) 
-H.c.) 



hi = cos{kx + ky) {a^dk + «l<^Lfc + H.c.) 
+ cos{kx - ky) {b^Ck + b^cij^ + H.c.) 



(12) 



(13) 



The first term in Eq. (ITOl ) is the zero-th order energy, denoted 
by £^0, which is independent of the angle (p at J = 1, as 
mentioned previously. By applying the Bogoliubov transfor- 



mation, we obtain a diagonalized form of the Hamiltonian as 

Hsw = Eo + AE+ Yl ^r'<'^"r'> (14) 

where a^ is the boson operator and subscripts r] and r]' take 
±. The energy dispersions are given as 



r/i2J,)=2SJx;i + r,'Y^ 



(15) 



with 



X^ = {J + 1 — (J — 1) cos 2ip} — 27] J cos kx cos ky cos^ Lp 
X {J + l-(J-l)cos2(/:)}, (16) 

and 

Y^ =2(cosV + -^sinV) 

X Y 4J^ sin^ kx s\v? ky cos^ Lp + (cos kx -\- r] cos ky)'^ sin^ (/:). 

(17) 

In the cases of(p = and 7r/2, the dispersion relations are 
reduced to 



^r7(2J.) = 2^^1 + wVcos(fe,+7//c^), 



(18) 



and 



.,^^ 



/(2 J^) = 2SJj {J + 7?'(cos /c^ + 7/ cos ky)/2}, (19) 



respectively. The second term in Eq. (O is a correction due 
to the zero-point vibration given by 



A^ = ^ 



- 2 J^N Sices'^ LP + Jsin^ (/:)). 



(20) 



In Fig. [71 the ground state energy including the quantum 
correction, Eq + A£^, is plotted as a function of J. Reduc- 
tions from the MF energies are remarkable around J = 1. At 
J = 1, energy for the limit of J ^ 1 — is lower than that 
for J ^ 1 + 0. Energies between the two configurations are 
shown in Fig. [8] as a function of the rotation angle Lp. En- 
ergy reduction due to the zero-point vibration is the largest 
at (/? = 0. That is, the (T/^ ^^) order is stabilized among 
the continuous degenerate configurations. This is attributed 
to the one-dimensional character of the spin- wave dispersion 
relation which gives rise to the large excitation density at low 
energies. 



B. Jordan- Wigner Method with MF Approximation 

In order to analyze the orbital state beyond the spin- wave 
approximation, we use a combined method of the Jordan- 
Wigner transformation and the MF approximation. We term 
this method MF+ID for simplicity. The MF approximation is 
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FIG. 7: Ground state energy at T = obtained by using sev- 
eral methods. Dotted and dashed-dotted lines represent results ob- 
tained by the MF approximation and the spin- wave method, respec- 
tively. Solid line represents the result by MD-i-lD method. Green 
broken line shows the coexistent area in the first-order phase transi- 
tion. Filled circles represent the results by the Lanczos method in a 
cluster of A^ = 32. 



given by 

i^MF-iD/(2 J.) = JY. ^^^/ -hY.T^+L {T^f , 

(21) 

where the Neel-type MF, (T^) = {TL^^), is assumed. A 
symbol < ij >' represents a NN ij pair in a one-dimensional 
chain, and L is a size of a chain. The transverse field is given 
by 



h 



-2 (T^) . 



(22) 



By introducing the Jordan- Wigner transformation, the PS op- 
erators are represented by a spin-less fermion operator Ci such 



asl^;^ 



c]ci 



\ and others. The fermion model is diagonal- 
ized by using the Bogoliubov transformation, and the follow- 
ing Hamiltonian is obtained; 



Hmf-ib/{2J^) = ^ ^fc ( aW -^]^L {T'Y , (23) 



where a^ is a fermion operator introduced by the Bogoliubov 
transformation defined by 
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FIG. 8: Ground state energy at J = 1 obtained by the spin- wave 
approximation (solid line). Broken line represents the MF energy. 
Schematic PS configurations assumed in the spin-wave approxima- 
tion are also shown. 



ak = UkCk - VkC 



with the coefficients given by 
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The eigen energy for the fermion is given as 



Eh 



J2 



Jhcosk + /i^, 



(24) 



(25) 



(26) 



where k is the wave vector along the (11) or (11) directions. 
We calculate the expectation value (T^) by solving the self- 
consistent equations in Eq. ([211 and Eq. (l22b . At zero temper- 
ature, we have 



{no = j^{cU)o-l 



L 

1 

2^ 



-h -\- ^ cosk 



Jh cos k -\- h"^ 



--dk, (27) 



applied to the NN interaction, T^Tj, and the one-dimensional 
chain under the MF is analyzed by the Jordan- Wigner trans- 
formation.^^ This treatment is justified for J <C 1. 

By introducing the MF approximation in the NN inter- 
action, the Hamiltonian for each diagonal chain is mapped 
onto the independent one-dimensional transverse Ising model 



where (• • • )o represents an expectation value at zero tempera- 
ture. The self-consistent equation is given as 
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FIG. 9: The left and right hand- sides of the self-consistent equation 
in Eq. (|28] ) at T = 0. Energy versus (T^) curves for each J are also 
shown. 



At finite temperature, the partition function and the free en- 
ergy are obtained from Eq. (1211 as 



^ = e-(2J.)L/?(T«>^-Q2cosh 



(2J.) 



/9-Efc 



(29) 
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FIG. 10: Phase diagram obtained by using the MF-i-lD method. A 
red circle represents the tricritical point. Broken line indicates the 
temperature below which the one-dimensional T^ correlation devel- 
oped. The lines in J < 1 and J > 1 are plotted in the different 
scales at left and right figure, respectively. A shaded area implies 
the coexistence region in the first-order phase transition. Energy ver- 
sus (T^) curves at (J,T/(2J^)) =(0.2,0.05), (1.1,0.05), (1.25,0.05) 
and (1.4,0.05) are also shown. The inset shows the extension around 
J= 1.2 and T/J^ = 0.1. 



and 



, 1J z Ej^ 

2 cosh — 

T 2 



(30) 



respectively. 

In Fig. [21 we plot the left- and right-hand sides of the self- 
consistent equation in Eq. ([28l) in a region of the positive effec- 
tive field (h > 0). When J is less than one, the equations have 
solutions at h = and at one positive value of h. At J = 1, 
slopes in the two curves coincide with each other at /i = 0. 
Additional solution appears in the case of 1 < J ^ 1.35. 
The energies are calculated as functions of (T^) (see the in- 
sets of Fig. [9]). The energy has three minima in the region of 
1 < J < 1.35. In particular, in 1.2 < J < 1.35, E shows 
absolute minima at (T^) = 0. A discontinuous change in the 
stable (T^) at J = 1.2 implies the first-order phase transi- 
tion. This is attributed to an existence of an inflection point 
in the h-{T^) curve, i.e. the magnetization curve in the trans- 
verse Ising model at T = 0. It is well known that this model 
shows a second-order phase transition at T = 0, a quantum 
critical point, at a certain value of h/J, where the magnetic 
susceptibility diverges due to an inflection point in the mag- 
netization curve. In this sense, the present first-order phase 
transition originates from competition between the directional 
dependent PS interactions, that is, TfTj along the NN bonds 
and TfTj^ along the NNN bonds in the compass model. 

Results at T = obtained by the ID+MF method to- 
gether with the results by other methods are summarized in 
Fig. [71 Broken line represents a region where solutions of 



(T^) = and (T^) ^ coexist, and an open circle indicates 
a point where the absolute minima change from (T^) 7^ to 
(T^) = 0. At J = 1, where the continuous degeneracy exists 
in the MF solutions, the (T^) order is realized in the ID+MF 
method. These results are consistent with the results in the 
spin- wave approximation. 

Next we present the phase diagram at finite temperature. 
We suppose an existence of the tricritical point at a certain 
(JjT), since the second-order phase transition appears at J = 
where the model is reduced to the two-dimensional Ising 
model, and the first-order phase transition is confirmed at T = 
as explained above. We expand the free energy in Eq. (l30l ) 
at the vicinity of (T^) = as 



F/{2J,L) = fo + {T'ff2 + {T'rh 
with coefficients given by 

/o = —tin 2 cosh 

/2 = 1 



1 



, ^ sech , ^ 
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4t 



(34) 



where t = T/{2Jz). The second-order phase-transition point 
is given by /2 = and /4 > 0, and the tricritical point is 
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FIG. 11: Amplitude of the {T^^^^)) order obtained by the MF+ID 
method (bold line) and that by the spin- wave approximation (dashed- 
dotted line). Broken line represents the results for the coexistent 
region in the first-order phase transition obtained by the MF+ID 
method. The square root of the orbital correlation function S^^ (tt, tt) 
calculated in the Lanezos method are plotted by filled circles. 



given by /2 = /4 = 0. The finite-T phase diagram is pre- 
sented in Fig. [TOl We also plot the first-order phase transi- 
tion points and the hysteresis region determined by the free 
energy in Eq. (l3Ql) . Broken line represents a crossover, be- 
low^ which the one-dimensional T^ correlation is developed, 
and is numerically determined by a peak in the specific heat 
C = -T{d'^F)/{dT'^) at (T^) = 0. This line does not 
depend on J, because the present model in Eq. (|2T1) is re- 
duced to the independent one-dimensional Ising model, when 
(T^) = 0. Both the results obtained by the MF+ID method 
and the results by the classical MC method (see Fig. [5]) show 
that the phase transition for {T^^ ,) is changed to be the sec- 
ond order to the first order through the tricritical point with 
increasing J. One discrepancy is seen at J ~ 1 in low tem- 
peratures; in the quantum phase diagram, the {TL ^) order is 
realized up to J ~ 1.2 even at T = 0. This is a kind of order 
by fluctuation phenomena due to the quantum fluctuation. 



C. Exact Diagonalization Method 

To examine the orbital state at T = in more detail, we 
adopt the exact diagonalization method based on the Lanezos 
algorithm. We use a 4>/2 x 4>/2-site cluster, where edges 
are parallel to the (11) and (11) directions, with the periodic 
boundary condition. Calculated energy is plotted as a func- 
tion of J in Fig. [71 The results are good agreement with the 
results obtained by the MF+ID method for J < 1, and with 
the results by the spin- wave approximation except for a region 
of J-1. 

The square root of the orbital correlation function 
S^^ (tt, tt) calculated by the Lanezos method is compared with 
the ordered moment of the (Tf. ^^) order obtained by other 
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FIG. 12: Temperature dependence of the binder ratio (see Eq. (|35}) 
obtained by the QMC method with L = 6 - 12 where 5 x 10^ MC 
steps are used for measurements. Statistical errors are estimated from 
64 independent runs. 



methods (see Fig.[TT]). The results obtained by three methods 
coincide with each other in a region of small J. However, at 
the vicinity of J = 1, large discrepancies between the three 
results are observed. Obtained ordered moment in the MF+ID 
method is larger than that in the spin- wave approximation. 
This tendency might be due to underestimation (overestima- 
tion) for the fluctuation in the MF+ID method (spin- wave ap- 
proximation). Data for the correlation function obtained by 
the Lanezos method are located between the results of (T^) 
by the MF-ID method and the spin wave approximation. 



D. Quantum Monte Carlo Simulation 

We present the numerical results obtained by the quantum 
Monte Carlo (QMC) method. Before showing the numeri- 
cal results, we touch signs of the exchange constants in the 
Hamiltonian in relation to the negative sign problem. Signs 
of the two exchange constants are positive in a view point 
of the perturbational calculation, but these signs can be re- 
versed by the following way. A checkerboard lattice is de- 
composed into the NN and NNN bond networks where T^ and 
T^ components are only concerned, respectively. Since the 
two networks are bipartite, signs of the exchange interactions 
can be reversed by introducing the unitary transformations of 
Tf -^ Uy^{7r)TlUy{7r) for the sites (ix.iy) of ix + iy=odd, 
and i;^ -^ U-^{Ti)T^Uy{Ti) for the sites of i^=odd. The 
unitary matrix t/^y (tt) represents the tt rotation around the T^ 
axis. In the simulations, we introduce the above transforma- 
tion, and the negative sign problem does not appear. In the 
manuscript, we choose signs of the exchange constants to be 
positive. 

We perform the continuous imaginary-time method with 
the loop algorithm in the ALPS library. ^^'^^ We use ^/2L x 
>/2L-site clusters (L = 6 — 16), where edges are parallel to 
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FIG. 13: (a) Detailed temperature dependence of the Binder ratio at 
J = 1.0 with L = 8 — 16. The arrows indicate the crossing points 
for the two curves in the L and L + 2 site clusters, (b) Scaling plot 
of the Binder ratio. We chose z/ = 1 and Tc/{2Jz) = 0.192. 




FIG. 14: Temperature dependence of the specific heat obtained by 
the QMC method, where 4 x 10^ MC steps are used for measure- 
ments. 



the (11) and (11) directions, with the periodic boundary con- 
dition. To calculate the physical quantities, 6x10^ — 5x10^ 
MC steps are used. Statistical errors are estimated from 4-64 
independent runs. 

In the region of J ^ 1, increasing of the correlation func- 
tion S^^ (tt, tt) at a certain temperature is observed (not shown 
in figure). The results indicate a possibility of the (T/^ ^x) or- 
der. In order to determine the critical point of this order, we 
utilize the Binder ratio defined by 






(35) 



In principle, this quantity does not depend on the cluster size 
at critical temperature. This is shown by utilizing the scaling 
relation given by 



g = fg\L'/''{T-T,) 



(36) 



where u is the critical exponent for the correlation length and 
fg is the scaling function. Figure [12] presents the temperature 
dependence of g for the several values of J and N. In the data 
sets for J = 0.25 and 0.75, the crossing points are observed. 
Deceasing of T at the crossing point with increasing J is con- 
sistent with the results obtained by the classical MC method 
and the MF+ID method (see Figs. [3l and [TOb. 

We focus on values of g at the crossing point, termed gc, in 
Fig. [121 It is known that, in general, gc does not depend on 
J, and about 0.85 for the two dimensional Ising universality 
classi^^ This figure shows that ^c's at J = 0.25 and 0.75 are 
close to this value, and the transitions are expected to belong 
to the two-dimensional Ising universality class. As for the 
case at J = 1, detailed results of ^ are presented in Fig.[T3la). 
Up to the results of L = 16, curves for different L do not cross 
with each other at same point. A value of g at the crossing 
point in the L and L-\-2 site clusters increases with increasing 



L, and might approach to 0.85 in the case of larger L. On the 
other hand, as shown in Fig. [T3lb), the Binder ratios plotted 
as functions of {T/Tc — 1)L^/^ are fitted by a single curve 
in the case of L > 12, and i/ = 1 and Tc/(2J^) = 0.192 
are obtained. The obtained value of u is consistent with the 
two dimensional Ising universality. From these analyses, we 
suppose that, at J = 1, the second order phase transition of 
{T^^ ^) is realized, and larger size clusters are required to 
examine gc than the clusters where the finite- size scaling for 
u works well. 

In the region of J ^ 1, on the other side, the calculated 
correlation function S'^^(7r, tt) does not show remarkable de- 
velopment with decreasing T. We examine the crossover tem- 
perature below which the one-dimensional T^ correlation de- 
velops for J ^ 1, as suggested in other calculation methods. 
The temperature dependences of the specific heat for several 
J and N are shown in Fig. [141 The maxima of the specific 
heats are indicated by small arrows. The crossover tempera- 
ture where C takes its maximum decreases with decreasing J. 
This tendency is similar to that observed by the classical MC 
method (see Fig. [5]), but is in contrast to that by the MF+ID 
method (see Fig. \T^ where the crossover temperature does 
not depended on J. This is attributed to the approximation in 
the MF+ID method in which the each one-dimensional chain 
is treated to be independent in this region. 

Next, we introduce careful examinations for the orbital 
states around J = 1. An accuracy of the MC simulation is 
checked by calculating the auto-correlation time for the auto- 
correlation function defined by 



C-(r) = {r(;,,)(r)T(;,,)(0)) 



{r(;,,)(r)){r(;,,)(o)>, 

(37) 



where Tf-^ ^■. (r) is the staggered orbital moment for the r-th 

configuration in the Markov chain. ^^ In the simulations where 
the MC steps are taken to be 10^, saturations for the auto- 
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FIG. 15: A scaling plot for the susceptibility at T/(2J^) = 0.15. 
We chose 7 = 7/4, zy = 1 and Jc/{2Jz) = 1.095. 
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FIG. 16: A scaling plot for the susceptibility at T/{2Jz) = 0.125. 
We chose 7 = 1, z/ = 5/9 and Jc/{2Jz) = 1.089. Inset shows 
a scaling plot where 7 = 7/4, v — 1 and Jc/(2Jz) — 1.088 are 
chosen. 



correlation times are observed above T/{2Jz) = 0.125, but 
not observed below T/(2 J^) = 0.1, and the numerical results 
obtained above T/{2Jz) = 0.125 are reliable. The orbital 
susceptibility at gr = (tt, tt) defined by 



X 



dr[{e 



tH 



%,nf 



-tH 



Tf^,,))-(Tf^,,))^], (38) 



is calculated at T/(2J^) = 0.125 and 0.15. The MC steps 
are chosen to be 1 x 10^ for T/(2J^) = 0.15 and 5 x 10^ 
for T/(2J2) = 0.125, and averaged values in the 4-times 
measurements are calculated. When temperature and system 
size are fixed, the susceptibility calculated as a function of J 
shows abrupt increase at a certain value of J which is termed 
Jc{T, L). When we assume that this points are the continuous 
critical points in the J— T plane, we expect a linear correspon- 
dence between J and T near the points. Therefore, from the 
conventional scaling form for the susceptibility as a function 
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FIG. 17: Phase diagram obtained by the QMC method. The lines 
in J < 1 and J > 1 are plotted in the different scales at left and 
right figure, respectively. Filled and open circles for the second-order 
phase transition are obtained by the finite size scalings of the Binder 
ratio and the susceptibility, respectively. Triangles for the crossover 
points are determined by the specific heat. 



of T, given by 



X- 



L^/-f^[L'/^T-T,)], 



(39) 



where 7 is the critical exponent and /^ is the scaling function, 
the following scaling relation as a function of J is expected 



X 



L^/-fJ[L^/-{J-J^ 



(40) 



where we introduce a scaling function f^. We suppose that, 
near the critical points, x/L^^^ versus L^^^{J — Jc) data for 
several L are on a single curve. 

A scaling plot at T/(2 J^) = 0.15 is shown in Fig. [T5l where 
(7, z/) = (7/4, 1), expected from the two-dimensional Ising 
universality class, and Jc = 1.095 are used. The optimized 
values obtained by the least-squares fit are (7, z^) = (1.5 ± 
0.8, 1.07±0.11),and Jc = 1.092±0.004. Scaling plot works 
well; numerical data obtained by several N are fitted by a 
scaling function. This analysis indicates that the second-order 
phase transition line continues from (T/(2 J^), J) = (0.28, 0) 
to (0.15, 1.092). 

On the contrary, the scaling analyses with the exponents 
(7, u) = (7/4, 1) do not fit the numerical data at T/(2 J^) = 
0. 125, as shown in the inset of Fig. [161 A different plot, where 
(7,z^) = (1,5/9) and Jc = 1.089 are used, is presented in 
Fig.[T6]for the data at T/(2J^) = 0.125. The optimized val- 
ues by the least-squares fit are (7, z/) = (0.8±0.5,0.58±0.08) 
and Jc = 1.094 ± 0.006. All data obtained in different N 
are almost fitted by a single function. The values (7, z^) = 
(1, 5/9) are the critical exponents for the two-dimensional tri- 
critical Ising universality class obtained by the c = 7/10 con- 
formal field theory j22i2S 

Phase diagram obtained by the finite- size scaling anal- 
yses in the Binder ratio and the susceptibility is given in 
Fig. [T71 As explained above, through the scaling analyses, 
we propose a possibility that the tricritical point exists around 
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(T/(2J^), J) = (0.125, 1.094). Because of an accuracy of 
the QMC simulation, the first-order phase transition expected 
below (T/(2J^), J) = (0.125, 1.094) is not confirmed by the 
numerical simulation. However, an existence of the tricritical 
point is reasonable by taking into account of the results ob- 
tained by other methods of the classical MC simulation and 
the MF+ID method shown in Figs. [3l and [TOl 



V. DYNAMICAL ORBITAL STATE 

In this section, we present numerical results for the excita- 
tion spectra in the checkerboard compass model. The excita- 
tion spectra are calculated by using the spin- wave approxima- 
tion, the MF+ID method and the continued fraction expan- 
sion method based on the Lanczos method. In the spin- wave 
approximation and the MF+ID method, the (Tf ^.) order is 
assumed. Results are presented in the Brillouin zone for the 
orbital disordered phase. 

The excitation spectra obtained by the spin- wave approxi- 
mation and the MD+ID method are explicitly given by 



UJ, 



(±) 

SW;i 



.j^/{2Jz) = V 1 - Jcos{kx ± ky), 



(41) 



and 



^MF: 



;J(2J.) 



J2 
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Jhcos{ka,±ky) ^h^, (42) 



from Eq. (fTSl) and Eq. (l26l) respectively. In the limit of 
J <C 1, ^sw-fc coincides with cjj^p.^ under the assumption 
of h = —2 (T^) = ±1. This is reasonable because the two 
approximations are equivalent with each other in this limit. 
Spin waves show one-dimensional character; dispersions ap- 
pear along (11) or (11) directions in the Brillouin zone. This 
is because, PS fluctuations in the {T? ^^ ) ordered state are 
caused by the interactions between T^ along the diagonal di- 
rections on the checkerboard lattice. In the Lanczos method, 
we calculate the dynamical correlation function given by 



S'\qM = 



-im(r^- 



1 



H^Ea 



IT] 






(43) 



where / = (x, ?/, z), £^^ is the ground- state energy and r] is 
an infinitesimal constant. In the numerical calculations, the 
system size is taken to be A^ = 32 and r] is chosen to be 
r]/{2J,) = 0.01. 

We show the excitation spectra for several J obtained by the 
three methods in Fig.[T8l where the momenta are varied along 
arrows shown in Fig. [19] Three results show good agreement 
with each other in the case of J = 0.5. Discrepancies between 
the three results are remarkable around J = 1. In particular, 
noticeable differences are observed in the lowest energy exci- 
tations; the lowest excitation energy by the spin wave approx- 
imation (the ID+MF method) is the lowest (highest) among 
the three results. A zero-energy peak in S^^{q) at gr = (tt, tt) 
is due to the static staggered correlation for T^ in the ground 
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FIG. 18: Spin wave dispersion relations obtained by the spin- wave 
approximation (dashed lines), and by the MF+ID method (dashed- 
dotted lines). Three dimensional plots for the dynamical PS correla- 
tion functions for the transverse component S^^{q, cj) + S^^{q, cj) 
(bold lines) and those for the longitudinal component S^^ (q, cj) (dot- 
ted lines) obtained by the Lanczos method are also shown. Param- 
eters are chosen to be (a) J = 0.5, (b)J = 0.9, (c)J =1.0 and 
(d)J = 1.1. 
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FIG. 19: The first Brillouin zone for the checkerboard lattice. Mo- 
menta in Fig. [18] are varied along the arrows. Bold lines represents 

kx^kz — ivr (see text). 



the Goldstone mode but is due to the linear spin wave approx- 
imation, and reflect the continuous degeneracy in the MF so- 
lutions at J = 1. We expect the dispersions are gapful when 
the higher order corrections in the spin wave approximation 
are taken into account. 

Let us focus on the results obtained by the Lanczos method. 
The present results in J = 0.5 well reproduce the results ob- 
tained by other two methods. With increasing J up to around 
J = 1, except for the lowest peaks, almost all peak intensi- 
ties are diminished and a number of small incoherent peaks 
appear. This might be attributed to the magnon-mangnon 
interaction which becomes remarkable when the system ap- 
proaches to J = 1. This result is related to the amplitude 
of the {T^^ ^^) order as well as the corresponding correlation 
function shown in Fig. [TTl where their reductions are due to 
the spin wave excitations. As for the lowest coherent peaks, 
even in the results by the Lanczos method beyond the spin 
wave approximation, their energies are almost flat along the 
lines of kx^ky = ±7r. These are shown in Fig.[T9]and corre- 
spond to the momenta for the PS configurations stabilized in 
J ^ 1.35. In this sense, the softening of the lowest coherent 
peaks implies a precursor of this PS configuration, although 
the phase transition at T = is of the first order. 



VI. DISCUSSION AND SUMMARY 

We discuss the present results in the checkerboard orbital 
compass model in comparison with the square lattice orbital 
compass model (SLCM). There are a number of theoretical 
studies in the orbital compass model on a square lattice de- 
fined by 



H 



J. Yl ^^^! ^J-Y. ^'^. 



<ij>z 



E 



(44) 



where the first and second terms are the NN interactions along 
the horizontal and vertical directions on a square lattice, re- 
spectively. There are the generators, which are similar to 
Eq. ([5]) in the present model, defined as Pi = ]\^^i Tf and 
Qm = Yli^m^i ' where / and m indicate the l-th row and the 



m-th column on a square lattice, respectively. The Hamilto- 
nian commutes with Pi and Qm for any / and m. The ground 
and excited states at T = have been studied by several meth- 
ods. It was shown in the anisotropic case, i.e. Jx ^ Jz that on 
a L X L-site lattice, the low energy spectrum consisting of 2^ 
states collapse exponentially fast with each other with increas- 
ing a system size. The first-order phase transition might occur 
at the symmetric point of Jx = Jz^ In contrast to SLCM, 
in the present checkerboard model, the two states realized in 
the large and small limits of J = Jx/{2Jz) are not symmet- 
rical with each other. The Neel-type symmetry-broken state 
is stabilized in the region of J ^ 1, and the 2^^-fold degen- 
erate staggered T^ ordered state along the diagonal directions 
appear from J = oo down to around J = 1. The first-order 
transition between the two occurs around J = 1.35 at T = 0. 
Just at J = 1, accidental continuous degeneracy is observed 
in the classical ground state in the present model as well as in 
SLCM. This is lifted by the quantum fluctuation and the Neel- 
type long-range order of T^ is realized at J = 1 and T = 0. 
In contrast to a number of studies at T = 0, little is known 
about the finite-T quantum states in SLCM. One of the reason 
is that any ordered phases are not expected to exist at finite 
temperature except for Jz = Jx-^^ There is an ordered phase 
in the present model, and the finite temperature phase diagram 
is obtained by the QMC simulation. 

By utilizing several methods, as well as QMC, we con- 
clude that there is a tricritical point around J = 1 at finite 
temperature. In the scheme of the ID+MF method, this is 
understood in analogy with the magnetization curve in the 
trans verse-Ising model, and originates from the directional 
depending interaction. The present results provide clue in- 
formation to reveal finite T quantum states in other-types of 
the compass models. It is also shown that, even on the ge- 
ometrical frustrated lattice, a conventional Ising model does 
not show a tricritical point. As an example, let us consider 
the Ising model on a checkerboard lattice where the inter- 
actions along the horizontal/vertical and diagonal directions 
are of JzS^Sj and JxS^Sj, respectively. There is a critical 
point for a certain value of Jx/ Jz^ termed Jc, where a macro- 
scopic number of degeneracy exists in the classical ground 
state. When the ID+MF method is applied to this model, a 
self-consistent equation corresponding to Eq. (|28l ) is obtained 
as a conventional MF type of (6'^) ~ tanh[(6'^)/T]. It is triv- 
ial that there is not a reflection point in this curve, in contrast 
to to Eq. ([28]) and a tricritical point is not expected. 

In summary, we study the orbital compass model on a 
checkerboard lattice where the interactions along the horizon- 
tal/vertical and diagonal directions are given as JzT^Tj and 
JxTfT^, respectively. The classical and quantum models are 
analyzed by several analytical and numerical methods. We 
obtain the finite temperature phase diagram as a function of a 
ratio of J = Jxji^Jz)- The Neel-type long-range ordr for T^ 
occurs for J ^ 1, and a crossover from a disordered state to 
the T^ correlated state along the diagonal chains is observed 
for J ^ 1. A reentrant feature of the Neel-type T^ order is 
shown around J = 1 due to the thermal order-by-fluctuation 
mechanism. A tricritical point around J = 1 is identified by 
QMC and ID+MF methods. This is understood from a view 
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point of the magnetization curve in the transverse Ising model 
and originates from the different types of the two competing 
interactions, i.e. JzTfTj and JxT^T^. Excitation dynamics 
in this model are also examined. In the vicinity of the phase 
boundary in the Neel-type T^ ordered state, the softening of 
the lowest coherent excitation peaks are confirmed along the 
lines of /ca^ ± /c^ = ±7r. This is interpreted as a precursor of 
the one-dimensional T^ order stabilized above J ~ 1.35. The 
present studies do not only provide new insights in a combi- 
nation of orbital frustration and geometrical frustration, but 
also help to reveal the finite T quantum states in other-types 



of orbital compass models. 
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